'''
Created on Sep 22, 2011

@author: oabalbin
'''

def read_positions(ifile):
    '''
    '''
    ifile = open(ifile)
    for l in ifile:
        if l.startswith('#'):
            continue
        fields=l.strip('\n').split(',')
    
    ifile.close()
    return fields

def create_rsID_list(fields):
    '''
    rs41356551|chrM@456|C>T
    '''
    rsIDs=set()
    for f in fields:
        rsIDs.add(f.split('|')[0])
    return rsIDs

def read_hapmap_rsID(file):
    '''
    '''
    ifile = open(file)
    rsIDs=set()
    for l in ifile: 
        rsIDs.add(l.strip('\n').split('\t')[0])
    ifile.close()
    return rsIDs

def main_rsIDs_comp(dbsnp,hapmap_rsd):
    '''
    '''
    fields = read_positions(dbsnp)
    rsIDs_dbSNP=create_rsID_list(fields)
    rsIDs_hapmap=read_hapmap_rsID(hapmap_rsd)
    
    rsIDisec=rsIDs_dbSNP.intersection(rsIDs_hapmap)
    print len(rsIDisec)
    
def main(file1,file2,file3):
    ifile = open(file1)
    ofile=open(file3,'w')
    ofile2=open(file3+'notfound_dbSNP','w')
    
    for l in ifile:
        if l.startswith('#'):
            continue
        fields=l.strip('\n')
    ifile.close()
    
    hfile =open(file2)
    for j in hfile:
        hf = j.strip('\n').split('\t')[0]
        freq=j.strip('\n').split('\t')[1]
        k=fields.find(hf)
        m=fields.find(',rs',k,len(fields))        
        feature = fields[k:m]
        if feature:
            rs,loc,change = feature.split('|')[0],feature.split('|')[1],feature.split('|')[2]
            chr,start, end = loc.split('@')[0],str(int(loc.split('@')[1])-1), str(int(loc.split('@')[1])+1)
            
            ol=[chr,start, end, rs+'|'+change+'|'+freq,'\n']
            ofile.write(",".join(ol).replace(',','\t'))
        else:
            ofile2.write(hf+'\n')
    
    hfile.close()
    ofile.close()
    
dbSNP='/exds/projects/alignment_indexes/gatk/hg19/dbSNP/dbSNP132_12_03_10_00-All_processed_04_05_11_3_db132.pos.base'
hapmap='/data/databases/hapmap/frequencies/freq/EUR/EUR_allele_freq_rsid.bed'
hapmap_out='/data/databases/hapmap/frequencies/freq/EUR/EUR_allele_freq_rsid_location.bed'
main(dbSNP,hapmap,hapmap_out)


